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ABSTRACT 



In this study a numerical model based on the hydrostatic Boussinesq equations is 
used to simulate atmospheric frontogenesis driven by an irrotational non-divergent 
deformation wind field. The equations are numerically integrated by using the semi- 
Lagrangian technique associated with two different time schemes: explicit and semi- 
implicit Both schemes produce realistic fronts after approximately 40 hours of model 
integration. Thesemi-Lagrangian semi-implicit scheme is more successful in handling 
the sharp gradients associated with the front. Also, the semi-Lagrangian semi-implicit 
equations are integrated with time steps as long as 3600 sec. producing solutions with 
relatively small errors. This indicates that this numerical scheme is appropriate for use 
in mesoscale regional models. 
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I. INTRODUCTION 



The development of fast computers has been of great importance for the 
atmospheric sciences. With them, scientists have been able to obtain a better 
understanding of the physical processes in the atmosphere and operational weather 
services have produced more accurate forecasts. But as pointed out by Robert (1981 ), 
scientists can also give their contribution by developing numerical algorithms that are 
computationally less expensive and at the same time accurate. This interaction between 
science and technology will allow the computational effort to be used in numerical 
models that will have higher resolution and that will be more accurate, as well, since 
more realistic and complex physical processes would be represented in these models. 

One of the atmospheric processes that has been under active study is the formation 
of fronts. In a study of non-geostrophic frontogenesis WLlliams(1972) reproduced 
numerically a realistic front, demonstrating the importance of nonlinear effects in this 
dynamical process. Since frontogenesis is a phenomenon that is characterized by the 
development of sharp gradients of the physical parameters it would be interesting to 
develop a numerical model based on techniques that can handle adequately the scale 
collapse problem. Kuo and Williams (1990) studied the use of different numerical 
techniques in the solution of a simple scale collapse problem. They concluded that the 
semi-Lagrangian method, developed by Sawyer (1963), would be adequate tor use in 
numerical models where strong gradients were present. 

The semi-Lagrangian technique has been used mostly in synoptic scale numerical 
models. Robert ( 1981) used this technique with the semi-implicit scheme (Robert et al.. 



1972) in a primitive equation model and obtained perfectly stable solutions using time 
steps 25 times larger than those that would be allowed by the Courant-Friedrichs-Lcwy 
<CFI ) stability criterion. Pudykiew'ic/ and Stamforth ( 19X4) examined the use of the 
scheme in several different conditions and concluded that besides its good stability 
properties, it is also accurate and flexible. 

The objective of the present study is to perform experiments using the semi- 
Lagrangian technique in a frontogenesis model and verify its suitability for use in the 
representation of atmospheric processes where sharp gradients are present. 

The semi-Lagrangian technique is introduced in the next chapter. The basic model 
equations are developed in chapter III. The numerical procedure employed for solving 
the equations is presented in chapter rV. In chapter V the experiments and results are 
discussed, followed by the summary and conclusions in chapter VI. 






II. THE SEMl-LAGRANGIAN SCHEME, 



This model uses a semi-L agrangtan algorithm similar to the one described by 
Robert ( 1981). 

The three time level semi-Lagrangian technique consists of evaluating the total 
derivative of a general dependent variable Q(r(f),t) following the trajectory of the 
fluid particle, using the approximation 

dQ _ Q(r(t»At),t+At)-<?(f(r-At),t-At) < 2 '> 

dt 2Ar 



where fit) represents the position vector of the fluid particle at time t. For sake of 
simplicity, constrain the problem to two dimensions, y and z In this case the position 
vector will be 

r(t) = y(t)f+z{t)k. (12) 

The location of the particle at the forecast time ft + A t) is chosen to be coincident 
with a grid point. Let the displacements of the particle along the y and z directions 
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during one time step be represented bv a and b . respectively If the grid point position 
is represented by 



f(r+Ar) = {y r z k ) 



the approximation (2.1) can be written as 

dQ _ Q(yj>z k ,t+ At)-Q(y } -2a,z k -2b,t- At) (2.4) 
dt 2A t 



The displacements a and b can be calculated iteratively using 

a n +1 - A t.v(y j -a H >z k -b n >t) , ( 2 5a ) 

b n * 1 = A t,w{y.-a n i z k -b n i t) , 

where the superscript a represents the o' iteration and v and tv represent the 
horizontal and vertical velocities, respectively. The initial estimates of a and b in the 
iterative process are defined by 
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a 0 = At.v(y.,z k ,t), 



(2 5c) 



b° - At.w(y r z k ,t) 



[n Pudykiewicz and Staniforth ( 1984) it can be found that the necessary condition for 
convergence of the iterative procedure described in (2.5) is 



At . max 



' dv 




dv 




dw 




dw 


\ 


k & 


y 


~dz 


9 


3y 


9 


dz 


/ 



< 1 



( 2 . 6 ) 



Furthermore. Kuo and Williams ( 1 990) point out that no more than three iterations are 
necessary for convergence to be attained. Robert (198 1) compared the use of two and 
four iterations in his primitive equation model and the results obtained were quite 
similar. 

Ln order to determine the values of the variable Q at time t-At it is necessary to 
use an interpolation scheme. In this model, bicubic spline interpolations (de Boor. 1 962) 
were used with the explicit formulations of the semi-Lagrangian scheme, using the 
algorithm described by Marchuck (1982). In the semi-implicit formulation only one- 
dimensional cubic splines in the y direction are necessary' to solve the problem: as will 
be shown in IV. C. 2. 
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III. BASIC EQUATIONS 



This model uses the hydrostatic Boussinesq equations, which assume 
incompressibility of the atmosphere. Friction and heating are neglected. It is also 
assumed that the atmosphere is bounded above by a rigid lid. Periodic boundary 
conditions are used in the horizontal. 

The Boussinesq equations can be written as 

* . AD 

dt dx 



dv 

~dt 




(3.2) 



^ + ^ + _ 0, (3.3) 

dx dy dz 



dd „ 
— - 0 , 

dt 



(3.4) 
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d<P _ g9 

dz % ’ 



where 



K - 




( 36 ) 



0 . r(-2.)‘-e 0 , 

Po 



(3.7) 



<P 



C p %(-?-) K + gz. 
Po 



(3.8) 



This model assumes that the Coriolis parameter /is constant, and that the bottom 
topography is flat. The variations in the x direction, along the front, are neglected, 
except for the <p field and the basic deformation field. 

With the assumptions above, equation (3.3) reduces to 

— + — = 0 (3.9) 

dy dz 



for departures from the deformation wind field. 
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With the use of a rigid lid assumption the pressure is not known at any level, and 
this will require some manipulation of the equations. 

Integrate the hydrostatic equation (3.5) from z = 0 to an arbitrary level z and 
obtain 



<P 




(3.10) 



The constant of integration, C, can be determined by taking the vertical average of 
equation (3. 10) and subtracting from equation (3.10) which gives 



z z 

<p = —fQdz- < — fddz> + <(p>, 
Vo 0 oo 



(3.11) 



where < > denotes the vertical average operator 



<F> 




(3.12) 



The constant of integration has been eliminated, but now an expression for <$>>must 
be found. 
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Apply the vertical average operator to the continuity equation (3.9) and use the 
vertical boundary conditions 

W = 0 at Z = 0 <3 13a) 

and 

w - 0 at z - H (3 13b) 



which yields 



d<V> 



= 0 . 



(3.14) 



Next, expand the total derivative in equation (3.2), multiply (3.9) by v and add 
to (3.2) to obtain 



dv dw dwv 

dt dy dz 




(3.15) 



Take the vertical average of (3.15), using the vertical boundary conditions (3.13), 
which gives 



a <v> a<w> a«p> v 

+ z_ -/<«> 

dt dy dy 
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Differentiating (3.16) with respect to v and using equation (3. 14) yields 



d 1 «p> _ _d 1 <w> fd<U> 
dy 2 = ay 2 



With appropriate boundary conditions, a solution for <4>>can be found. 

[n order to solve equation (3.1), an expression for — — must be obtained. 

dx 

Differentiate (3. 1 1 ) with respect to x and use the assumption that £is not a function of 
x which gives 



8<p _ d< <p > 
dx dx 



An expression for > must be obtained. Take the domain average of (3.2) 

dx 



{<^>} - 

dt dy 



(3.19) 



where {( )} denotes the horizontal average operator 




2 



(3-20) 
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Expand the first term of ('3.19) and use the vertical boundary conditions (3. 13) which 
yields 



{< 




$< v>} 

dt 



d<v> 

dt 



(3 21) 



Using periodic boundary conditions in they direction, it can be shown that the 
second term in (3.19) is zero. With the assumption of constant f, and using (3.21), 
equation (3.19) becomes 



— — - ~f[<U>). (3-22) 

dt 



Likewise, take the domain average of (3.1) and use the vertical and horizontal boundary 
conditions which yields 

cH<K>} cH<<p>) f (3 23) 

— — +/<v>. 

dt dx 



Equations (3.22) and (3.23) show that if {<u>} is set initially equal 

<v> 



will remain constant in time and so w ill ^ . Furthermore, if 

.{< > , dx 

zero initially, — ^ — will be identically zero for all time. 



dx 



to zero, <v> 
is also set to 
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In this study it will be further assumed that 
direction (y), which allows to write 



d<<p> 

ax 



is constant m the cross frontal 



di<<p>} d<tp> 

dx dx 



(3 24 ) 



With the assumptions above the following relation is obtained 

dtp _ d<(p> _ q 

dx dx 



for initial conditions 

(<«>} = <V> = 0. ( 326 > 



The dependent variables will be expressed as a combination of a basic and a 
disturbance part as follows: 



u(x,y,z,t) - U d (x,y) + u'(y,z,t ) , 



( 3 . 27 ) 



v(x,y,z,0 - F rf (jx,>) + v / (y,z,f), 



( 3 . 28 ) 



w(y,z,0 - w'( y>z,t). 



( 3 . 29 ) 
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90^,0 = 9 



<1 3(1) 






where the quantities with prunes represent the disturbance part of the variables 

Ihe basic wind field that drives the frontogenesis is the nondivergent and 
irrotational deformation field given by 

D, 

U d (x,y) sinh(fix)sin(|iy) , (3 32 > 

D, 

V d (x>y) = - — cosh (fijc) cos (jiy) . (3 33) 

r his w ind deformation field was originally used by Stone (1966) in his study of quasi- 
geostrophic frontogenesis. 

It has been assumed that the perturbations are independent of x If the 

expressions for u and v (3.27 and 3.28 respectively) are substituted into the basic 

equations, — and — will be functions of x. In this study, where the mam interest is 
dt dt 

in ex aminin g cross front variations, the equations will be evaluated at x - 0. This w ill 
make the mathematical formulations compatible w ith the previous assumptions and it 
is expected that the results will be physically consistent. 
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Therefore, the basic wind deformation field will become 



V d ( 0,y) - 0, 



(3.34) 



^(0,y) - -i^cosny. 



(3.35) 



The basic geopotential field # is chosen such that it will be in geostrophic balance 
with the wind deformation field as follows: 




(3.36) 




(3.37) 



Substitution of (3.36) and (3.37) into (3.27) and (3.28) respectively, and use of the 
expressions obtained for the dependent variables in equations (3.1)-(3.4) yields the 
following set: 



du 

dx 



-/v, 



(3.38) 
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(3 39 ) 



dv 

~dt 



- 4^ -/« - (V + v d ) 

dy 



dy 



J 



*♦^- 0 . ( 3 . 40 ) 

dy dz 



de 

dt 



0 . 



( 3 . 41 ) 



where the prunes were dropped from the disturbance variables. 

Finally, an easier way to obtain will be described. The dependent variable q> 

dy 

can be expressed as 



<P = <P 0 +< P,> 



( 3 . 42 ) 



where 



Vs 




( 3 . 43 ) 
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Differentiate (3.42) with respect to y to obtain 



3<p 


- * f dd dz. 


( 3 . 44 ) 


dy 


dy 0 O ^ dy 





The value of q> 0 is not known, by virtue of the rigid lid assumption In practice v 

<9<p n 

will be predicted using equation (3.44) with — 2 assumed initially equal to zero. After 

dy 

that, the v field will be adjusted in order to satisfy < v > = 0. This procedure will be 
equivalent to solving equations (3.11) and (3.17) for q > , differentiating the values 
obtained with respect to y, and substituting into (3.39). Equations (3.38)-(3.4l) and 
13.44) give a complete set that can be integrated in time in order to predict the values of 
the dependent variables. 
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IV. NUMERICAL SOLUTION 
A. INITIAL CONDITIONS 

The initial conditions used in this model are similar to the ones used by Williams 
et al.(1991) in their study of effects of topography on fronts. 

The initial temperature field is given by 

<30 if 

Q(y>z,0) - — - e (4 -*J 

oz 2 



where 



e 



- A cos(tiy)cos? 



[M 


n 


.UJ 


7 



zzz t ( 4 - 2 ) 



e - 0, z>z { (4 *) 

dd 

In the expression above — i is a constant. A is the amplitude of the temperature 

dz 

disturbance and z t is the height of the upper limit where the temperature disturbance is 
present. This initial temperature field will confine the frontogenesis to the lower layers 
of the atmosphere, which will make the results more realistic. 
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The initial u component field is chosen to be in thermal wind balance with the 
temperature field and is given by 



u(y,z, 0) 



- s in(n_y) . X 

2/e 0 



(4.4) 



where 



X - (z-z,)- — sin 

1C 



Z-Zf) 



[\ 



it 



zzz. 



(4.5) 



X - 0, 



z>z t 



(4.6) 



and the condition 



«(y,z r >0) - o 



(4.7) 



has been used. 
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The initial vand tv fields are obtained from the solution of the quast-geostrophic 
circulation equation (Williams, 1972) 

2g_ 39 (4 

f 2 6 0 * dy 2 dz 2 / 2 0 o 3y dy ’ 

where 

V - - ii (4.9, 

dz 



w 



it 

dy 



(4,10) 



The vertical boundary conditions for the solution of (4.8) are set to 

i|/(z - 0) - i|f(z - H) - 0 (4.1 1 ) 



It can be shown, from the domain averaging of equation (4.2) that {< u>} will 
be zero initially. Also, the condition (4. 1 1) will allow (3.25) and (3.26) to be satisfied in 
order for the numerical solution to be consistent with the assumptions of the model. 
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B. THE SPATIAL GRID 

The numerical calculations are performed on a staggered grid, both in the 
horizontal and in the vertical directions. In the horizontal the variables are staggered 
following the scheme B described in Arakawa and Lamb (1977). This scheme, as shown 
in Haltiner and Williams (1980), allows a better representation of geostrophic 
adjustment processes. The variables are also staggered in the vertical. Pielke ( 1984) 
points out that the staggering of the u and w components in the vertical gives better 
solutions for the vertical velocity from the continuity equation. 

The horizontal grid uses I grid points separated by a constant grid space Ay. The 
vertical grid has K grid points and also uses a constant grid spacing, referred to as Az. 
Figure 1 shows the arrangement of the variables on the grid. 

C. FORMULATION OF THE MODELS 

1. Semi- Lagrangian Explicit Scheme 

In the semi-Lagrangian explicit formulation of the present model, the total 
derivatives will be approximated using the technique described in chapter II. The 
remaining terms of the equations are evaluated at the positions of the fluid particles at 
time level t. Introducing those approximations into the prognostic equations (3.38), 
(3.39) and (3.41) the following set is obtained, after dropping subscripts / and j . 

2At 



f v(y-a,z-b,t), (4 12 ) 
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• W 




• W 




• w 


z ( k + 1 ) 


• 0,<p 


• 

U,V, v d 


• 0,<p 


• 

u,v,V d 


• 0,<p 


z(k+£ ) 


• W 




• w 




• w 


z(k) 


• 0,<p 


• 

u,v,V d 


• e,q> 


• 

u » v *V d 


• 0,q> 


z{ k-j ) 


• w 




• w 




• w 


z ( k— 1 ) 


y(j-i) 


y(j-£) 


y<j) 


y(j+l) 


y(j+i) 





Figure 1: Staggered spatial grid 
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v (y>z,f + A 0 - v(y-2a,z -26,r- Ar) 

2At 




(y-a,z-b,t)-f.u(y-a,z-bj) 



-[v{y-a,z-b,t) + V d (y-a)] 



d l* 

dy 



(y-a), 



(4.13) 



9(y > z > r+Af)-Q(y-2^ > z-2^,f-Ar) _ Q 
2A t 



(4.14) 



The values of w&nd 47 are obtained from the diagnostic equations (3.40) and (3.45), 
respectively. 

2. Semi-Lagrangian Semi -Implicit Scheme 

The semi-Lagrangian semi-implicit scheme was introduced in this model 
following the development described by Robert et al.(1985). The scheme consists of 
separating the temperature into a basic state, dependent only on the vertical coordinate 
z, and a disturbance part. Next, an implicit treatment is given to those terms related to 
gravity waves, namely: horizontal pressure gradient, divergence and vertical motion. 

To obtain the system of equations in the semi-implicit formulation take the 
advective terms out of the total derivatives in equations (3.38), (3.39) and (3.44); 
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(4 15) 



d u u du 

+ H’ 

di dz 



- /v. 



V dv 

+w 

dt dz 



8<p 

dy 



dV, 

-f u ~( v +V d )—^- , 

dy 



(4.16) 



dt 



ae n 

+ w — - 0, 
dz 



(4.17) 



where the terms with subscript H represent horizontal components of the total 
derivative. 

Replace the temperature and geopotential fields by basic state and 
disturbance parts as follows: 

0(y,z,r) - 0*(z) + 0 / (y,z,f), (418) 



<p (y>z,t) - q>*(z) + <p / (y.z»*)» 



(4.19) 



where the superscripts (*) represent the basic state parts of the variables. 

Now introduce (4.18) and (4. 19) into equations (3.5), (3.40) and (4. 15)-(4. 17), 
and take the time average of those terms related to gravity waves, which yields 
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d H U 

dt 



+ w- 



du 

dz 




d H V dv 

+ w 

dt dz 



<v 

dy 



-fu-(v+V d ) 



3K, 

dy 



9 



dv dw 
dy dz 



d $ ae' 

i -W 

dt dz 




- 0, 



W m gV_ 
dz 0 ° ’ 

dz " 0 O * 



where ( ) denotes the time average operator. 



(4.20) 



(4.21) 



(4.22) 



(4.23) 



(4.24a) 



(4.24b) 
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Define the explicit terms at time t as 



du , 

r, 



<4 25) 



5v SV, 

r, - w— +/u + (v+KJ — - 
2 dz * dy 



(4 26) 



r 3 = w 



d&_ 

dz 



(4.27) 



Using (4.25)-(4.27) in (4.20M4.23), the following set of equations is obtained: 



^H U 

dt 



+ r. 



0 , 



(4.28) 



dt 



dtp 

dy 



+ r„ 



- 0 , 



(4.29) 
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(4.30) 



dv + dw 
dy dz 



— 50 * 

+ w 

dt dz 



+ r. 



= 0 . 



(4.31) 



Next, define the auxiliary variables 

F* - F(y,z,r+ At), ,4 - 2 “> 



F° - F(y-a,z,t), 



(4.32b) 



F~ = FO?-2<z,z,f- Ar), 



(4.32c) 



where 



a n +1 - Af.(v 0 (y-a'\f) + F </ °(} , -tfV))- ( 4 - 33 ) 



Note that the use of the time average of the vertical velocity w eliminates the 
need of considering vertical displacements in the calculation of the trajectories of the 
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fluid particles. Therefore, the interpolation will be performed only in the horizontal 
direction using one dimensional cubic spline functions. 



Define the following approximations for the time derivatives 
average respectively 



d H F _ F* -F~ 
dt 2At ’ 



r F* + F 
F = 



Use of approximations (4.34) and (4.35) into (4.28)-(4.3 1) gives 



U + ~U 

2 At 



+ r x = 0, 



V - V 

2 At 



1 

+ — 

2 



/ 0<p" + dtp ^ 






dy dy 



+ r ? - 0, 



and time 



(4.34) 



(4.35) 



(4.36) 



(4 37) 
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dv* 

dy 



dw 
+ 

dz 



0 , 



(4.38a) 



dv dw 
+ 

dy dz 



(4.38b) 



0 + -0* 1 00*/ + n / 4 , 9 , 

+ — iw + w ) + r a - 0 . (4.4y) 

2Af 2 dz 3 

Group the terms at time (t-At) and define new auxiliary variables p t ; 

( 4 - 40 ) 

p 2 - - v~ + Ar — - , (4 - 41 ^ 

dy 

p - + ( 4 4 2 ) 

dz 
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Similarly, group the terms at time ft t J/Mnd define the auxiliary variables 






<h 




(4 43) 



q 2 = v + + At 




(4.44) 



*3 



rd+ A * 

(j + At w 

dz 



(4.45) 



Use of (4.40)-(4.45) allows to rewrite (4.36), (4.37) and (4.39) as follows: 

q x - -p { -2 At r, , (4.46) 

q 2 = -p 2 -2At r 2 , (4.47) 

q 3 - -p 3 -2At r 3 . (4.48) 



The variables p, and r i can be calculated explicitly and with their values the 
q, variables can be obtained. 
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To calculate the values of the variables u. v. xv. 0 and <p at the time level 
[t t At) it is necessary to solve the system of equations composed by (4.38b), (4.43)- 
(4.45) and the hydrostatic equation (4.24b) The procedure used by Robert et al.( 1985) 
consists of first solving an elliptic equation for <p. 

fn this model it is not possible to state exact boundary conditions for <p . due 
to the use of the rigid lid assumption. However, the w field has the exact boundary 
conditions (3. 13). The following elliptic equation for wean be derived from the system 
of equations at time level ft + At ) : 

d 2 ^* + A t\g ( 36 * _ A tg #<h + ^2 ( 4 49) 

dz 2 & ty 2 c^y 2 dydz 

Once w was calculated, 0 was determined using (4.45), and - and v were 

d y 

obtained using the procedure described in Chapter III. After adjusting the v field to 

satisfy < v> = 0, a new field was calculated using the new values of v in (4.44) 

dy 

and the w field was also adjusted to satisfy the continuity equation (4.38a). 

Finally, both in the explicit and the semi-implicit schemes the spatial 
derivatives are approximated by centered in space finite differences. Following Monk 
( 1989), the derivatives are first evaluated at the grid points and, after that the values are 
interpolated for the position of the fluid particle. This procedure is computationally 
more efficient than the calculation of the derivatives at the fluid particle points, which 
would double the number of interpolations. 
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V EXPERIMENTS AND RESULTS 



A. PARAMETERS USED IN THE EXPERIMENTS 

I he iollowing values for the constants were used in all experiments 
L - 3600 km. 

H = 12 km. 

A = 12° K, 

0 = 300° K, 
g -• 9.8 lms‘, 
f = 10 4 s' 1 , 

D 4 = 10 5 s 1 . 
z t - 9 km, 

—i = 4° K km ' . 

dz 

For the semi-Lagrangian semi-implicit (SLSI) case there were two control runs 
using A y - 20 km and A z - 167 m. All the remaining experiments used normal 
resolution, with A y - 40 km and A z- 333m. The cross-sections of the initial fields are 
shown in Figs. 2-5. There is a warm front with thermally indirect circulation near v- 
900 km and a cold front with termally direct circulation near v- 2700 km. 

In order to assess the evolution of the frontogenesis process, the following 
parameter is defined at the lowest computational level. 



d - 



'A 61 

ae 



ay 



max 



(51) 
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Figure 2. Initial u component field (m/s). 




Figure 3. Initial k component field (m/s). 
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Figure 4. Initial ^component field (m/s). 




Figure 5. Initial 0 field (°K). 
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where A0 is the maximum variation of 0 along the y direction. This parameter, used 
by Williams et al.(199l), gives a reasonable measure of the width of the frontal zone 

During the first experiment wit' he semi-implicit scheme, a high frequency 
oscillation was observed in the frontal width. That noise was not present in the 
experiments carried out with the semi-Lagrangian explicit (SLEX) scheme. In order to 
eliminate the noise, a time filtering technique similar to the one described by Asselm 
( 1972) was used. The time filtering of a variable Q is defined as 

W) - <2(0 + YlQ(f + Ai)-2<?(r)+<3(t-Ar)] , ,52 > 

w'here Q(t+A /)has been previously obtained by using the respective predictive equation 
w ith the unaveraged value of Q(t). The time filter defined above is expected to affect the 
higher frequencies only. 

Several experiments were performed with the objective of assessing how the use 
of time filtering would affect the solutions, as will be described next. 

B. DESCRIPTION OF THE EXPERIMENTS 

The experiments basically consisted of carrying out time integrations of the model 
using the SLEX and SLSI schemes, starting from the same initial conditions and using 
different combinations of time step At and filtering factor y. The experiments were 
intended to give information about the following points. 

• How well the different schemes handle the formation of large gradients; 

• How the use of time filtering affects the final solutions; and 

• How the solutions are affected by the use of longer time steps. 
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Table I shows the experiments performed with the SLHX scheme, as a function 
of time step A t and filtering factor y. The experiments carried out using the SI. SI 
scheme are identified on Table 2. 

TABLE 1 EXPERIMENTS WITH THE SLEX SCHEME. 



A t (sec.) 


y = 0.0 


Y = 0.01 


y = 0.03 


180 


SLEX01 


SLEX02 


SLEX03 


360 


SLEX04 


SLEX05 


SLEX06 


600 


SLEX07 


— 


SLEX08 



TABLE 2. EXPERIMENTS WITH THE SLSI SCHEME. 



| At (sec.) 


-< 

11 

o 

o 


y=0.01 


Y=0.03 


y=0.05 


y=0.07 


y=0. 1 


I 180 


SLSI01 


SLSI02 


SLSI03 


SLSI23 


— 


- 


I 360 


SLSI04 


SLSI05 


SLSI06 


SLSI24 


— 


.... 


| 600 


SLSI07 


SLSI08 


SLSI09 


SLSI25 


— 




1 1200 


— 


SLSI10 


SLSI11 


SLSI12 


— 


.... 


1800 


— 


SLSI 13 


SLSI 14 


SLSI15 


— 


.... 


I 2400 


— 


SLSI 16 


SLSI 1 7 


SLSI18 


SLSI 19 




j 3600 


— 


SLSI20 


SLSI21 


SLSI22 


SLSI26 


SLSI27 
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There were two high-resolution experiments, both using the SLSI scheme, with 
time steps of 180 secs. In the first of them , identified by SLHROI, the time filtering 
factory = 0.01 was used . and in the second one(SLHR02) y was set equal to 0 05. The 
velocity components and temperature fields at time t - 40 hours, for experiment 
SLHROI are shown in Figs. 6-9 The occurrence of frontogenesis for the surface cold 
front and frontolysts for the warm front can be seen. These processes are due to the 
relation between the wind deformation and temperature disturbance fields, that gives 
convergence in the cold front region and divergence for the warm front. The wand v 
fields give a direct circulation around the cold frontal zone and an indirect circulation 
about the warm frontal zone. The u field develops cyclonic shear at the cold front. 

The evolution of the frontal width ( d-value ) for SLHROI is presented in Fig. 10. 
An almost linear decrease in the frontal scale can be seen. This curve is different from 
that obtained by Williams et al.(1991), because in that case there was a momentum 
diffusion term that gave an asymptotic behavior in the evolution of the frontal width 
towards a balance condition. In the present model a continuous reduction of the frontal 
scale is expected until a limit, when the grid resolution is reached. Figure 1 1 shows the 
evolution of the d-value for experiment SLEX0 1 . It can be seen that the minimum value 
of the frontal width is reached at approximately t - 40 hours. Note that this curve is 
more linear than that obtained in experiment SLHROI. Also, the d-value curves 
obtained for experiments SLEX02 and SLEX03 were virtually coincident with the one 
for SLEX01 . This result shows that the use of stronger filtering effect (larger y ) does not 
affect the evolution of the frontogenesis process represented by the SLEX scheme. 

During the first experiment with the SLSI scheme (SLSI01) the d-value curve, 
shown in Fig. 12, presented a high frequency oscillation. The time integration could be 
carried out until time t— 36 hours and after that, the iterative procedures no longer 
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y(km) 



Figure 6. u component field (m/s) for SLHR01 at t = 40 hours. 




Figure 7. v component field (m/s) for SLHR01 at t - 40 hours. 
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Figure 9. 6 field (°K) for SLHR01 at t = 40 hours. 
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Figure 10. Evolution of the frontal width (d-value) for SLHR01. 




Figure 11. Evolution of the frontal width (d-value) for SLEX01. 
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achieved convergence and the integration was terminated. It is important to note that 
the irregular behavior of the curve was not related to continuous amplification of the 
solution due to numerical instability and that the curve presents a general tendency to 
the frontogenesis solution obtained in the SLEX experiments. The high resolution 
experiment SLHR01 also shows a low frequency oscillation about its linear decrease in 
d. 



o 




Figure 12. Evolution of the frontal width (d-value) for SLSI01. 



In order to eliminate the noise the time filter described in the last section was 
introduced. Figure 13 shows the evolution of the frontal width obtained in experiment 
SLSI02, with the filtering factor y set equal to 0.01. It can be seen from the figure that 
although the filtering effect was small, it was adequate to suppress the high frequency 
oscillations. 
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Figure 13. Evolution of the frontal width (d-value) for experiment SLSI02. 



Figures 1 1 and 13 show that both for the the SLSI and SLEX experiments the 
minimum frontal width is reached approximately at time /= 40 hours and after that the 
d-value curves oscillate. Those oscillations can be related to the fact that the front 
reached a width corresponding to the model resolution and after that point the solutions 
were no longer physically consistent. 

Tables 3 and 4 show the minim um value of the frontal width for the SLEX 
experiments with time steps 180, 360 and 600 sec., and the corresponding SLSI 
experiments for which the time integrations could be carried out until t - 40 hours. 
Also presented is the time when the minimum d-value ocurred. 

Comparison of the values presented in Tables 3 and 4 show that the minima d- 
values obtained in the SLSI experiments are smaller than those obtained for the 
corresponding SLEX experiments with same A /and y. Also, in the SLEX experiments 
the minimum frontal width increases as A /increase whereas in the SLSI case there is no 
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significant change in the minimum d -value due to changes in the time step It is 
interesting to note that for the SLEX case a shorter time integration is necessary for the 
minimum to occur for increasing tune step, while the inverse tendency is observ ed in the 
SLSI case that is. a longer time integration is necessary for the minimum frontal width 
to occur, for larger A t. 

TABLE 3. MINIMUM D-VALUE FOR SLEX EXPERIMENTS. 



A t (sec. ) 


r 


Experiment 


Minimum 
d-value (km) 


Time (hour) 


180 


0.0 


SLEX01 


49.0 


45.5 


180 


0.01 


SLEX02 


49.0 


45.5 


180 


0.03 


SLEX03 


49.0 


45.5 


360 


0.0 


SLEX04 


68.4 


39.6 


360 


0.01 


SLEX05 


68.3 


39.9 


360 


0.03 


SLEX06 


68.3 


39.6 


600 


0.0 


SLEX07 


69.9 


39.3 


600 


0.03 


SLEX08 


69.8 


39.5 



Therefore, it can be concluded that the SLSI scheme gives a better representation 
of the scale collapse process than that obtained from the SLEX scheme, since narrower 
frontal widths can be resolved by the former than by the latter. 
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TABLE 4. MINIMUM D-VALUE FOR SLSI EXPERIMENTS 



A t (sec.) 


Y 


Experiment 


Minmum 


Time 








d-value (km) 


(hours) 1 


180 


0.01 


SLSI02 


43.4 


41.5 


180 


0.03 


SLSI03 


43.4 


41.7 


360 


0.01 


SLSI05 


43.4 


41 5 


360 


0.03 


SLSI06 


43.4 


42.0 


600 


0.03 


SLSI09 


43.1 


42.0 



Experiment SLSI03 used y= 0.03 and the d-value curve obtained was practically 
coincident with the one obtained in SLSI02, with y=0.01. 

In experiments SLSI04 and SLSI07 no filtering was used ( y =0.0) and 
approximately at time t- 37 hours the numerical integration was interrupted because 
the solutions were no longer able to achieve convergence in the iterative procedures 
However, the use of larger values of y (SLSI05. SLSI06. SLSI24, SLSI08. SLSI09. 
SLSI25) allowed the normal execution of the time integrations. Figure 14 shows the d- 
value curves for experiments SLSI08 and SLSI09. It can be seen that the oscillations 
were not completely filtered out with y =0.01 and a stronger filtering effect (y=0.03) was 
required. The experiments with longer time steps showed that as the time step increased 
a larger value for y was required for an appropriate filtering of the noise. Also, all 
experiments with the SLSI scheme showed that the largest amplitudes of the noise 
occurred during the first 24 hours of integration. Furthermore, those experiments 
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showed that the noise was not related to unbalance in the initial conditions, since the 
mass and wind fields were initially adjusted using the quasi-geostrophic circulation 
equation (4.8) and the oscillations were not present before t- 02 hours. The presence 
of the noise was investigated by examining the u, v and 0 fields obtained at the lowest 
computational level of experiment SLSI07. At time t =13 hours, that corresponds 
approximately to maximum amplitude of the noise, there was no evidence of small scale 
spatial oscillations in those fields. Although the origin of the noise could not be 
determined, the experiments showed that the use of the time filter was effective in 
e limin ating the oscillations with relatively small values of y. It is expected that this 



filtering would not significantly affect the lower frequency solutions. 




Figure 14. Evolution of the fronal width (d-value) for experiments SLSI08 and SLSI09. 



Another point of interest in this study was the computational effort employed in 
each of the experiments, since one of the objectives of the SLSI scheme is to allow stable 
solutions with long time steps. In order to evaluate the relative computational 
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efficiency, the CPU tune spent tor running each experiment for 42 hours was measured 
and normalized by the CPU time used in SLSI01. It was observed that the CPU time 
was not affected by changes in the value of y . so that the results are presented in I able 
5asa function of numerical schemeand timesteponly. As expected, the computational 
effort decreases almost linearly for increasing time step Also, for a given time step, the 
SLEX scheme is more efficient than the SLSI one. This is observed because the SI SI 
technique requires the solution of an elliptic equation every time step. It should also be 
pointed out that in this study the SLSI formulation makes the spatial interpolations 
necessary only in the horizontal direction, whereas in the SLEX scheme the 
interpolations are performed using bicubic splines. Thus, a formulation of the SLSI 
scheme using two-dimensional interpolations could give a further increase in the 
computational effort. The advantage of the SLSI scheme appears clearly when 
integrations with time steps as long as 3600 sec. are performed with perfectly stable 
solutions, while the SLEX scheme did not allow time steps longer than 600 seconds 
In order to assess the accuracy of the several solutions, a second high resolution 
(SLHR02) control run was carried out with the SLSI scheme. The parameters used 
were A t - 180 sec. and y = 0.05. This value for y was chosen to guarantee that the 
solutions would not be affected by noise. The solutions of the control run were linearly 
interpolated for the normal resolution grid point positions when necessary, because of 
the staggering of the variables. The differences between the values obtained for the 
control run and those obtained for the experiments with normal resolution were 
calculated at each grid point. The results were used to calculate the domain RMS 
differences for each variable. Two sets of comparisons were performed: In the first one 
the velocity components and temperature fields for selected SLSI experiments at t- 36 
hours were compared with the ones obtained from SLHR02 at the same time. That 



45 



time was considered appropriate to give a well characterized atmospheric front and 
sufficiently away from the inconsistent solutions obtained after the model resolution 
was achieved. 



TABLE 5. NORMALIZED CPUTIME FOR 42-HOUR MODEL INTEGRATION. 



| A t (sec.) 


SLEX 


SLSI 


I 180 


0.84 


1.00 


I 360 


0.42 


0.50 


1 600 


0.25 


0.30 


1200 


— 


0.16 


1800 


— 


0.11 


2400 


— 


0.08 


3600 


— 


0.06 



The experiments chosen for this first set of comparisons were those with y= 0.05 
for all time steps, except A t= 3600 sec., where y- 0.07 was used. The accuracy of the 
solutions was expressed in terms of RMS differences in the values the of u, v.w&nd 0, 
given in Table 6. 

The values obtained show that the u and 0 fields have small errors as compared 
to the range of the values observed in the domain (-35 to +22 <> K for 0 and -50 to + 15 
m s ' for u). On the other hand, the errors obtained in the v and w fields were relatively 
large compared with the magnitude of the reference values obtained for SLHR02. The 
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large differences can be associated with the fact that the r and » components are 
closely related to the divergence in the model and their values can be significantly 
affected by the presence of gravity waves. Time series of v and » were examined and 
they confirmed the presence of internal gravity waves in the solutions llius, a single 
sample would not be sufficiently representative of the actual level of accuracy obtained 
m those fields. 

TABLE 6 RMS DIFFERENCES FOR TIME T= 36 HOURS. 



Experiment 


At (sec.) 


Y 


u(m/s) 


v(m/s) 


w(nvs) 

*10' } 


0(°K) 


SLSI23 


180 


0.05 


0.170 


0.134 


0.58 


0.087 


SLSI24 


360 


0.05 


0.180 


0.168 


0.92 


0.088 


SLSI25 


600 


0.05 


0.220 


0.256 


1.62 


0.098 


SLSI12 


1200 


0.05 


0.373 


0.526 


3.40 


0.139 


SLSI15 


1800 


0.05 


0.554 


0.796 


5.87 


0.191 


SLSI18 


2400 


0.05 


0.764 


1.038 


7.23 


0.254 


SLSL22 


3600 


0.07 


1.204 


1.270 


8.56 


0.377 



The differences in the u and 0 fields obtained at each grid point for experiments 
SLSI23 and SLSI15are shown in Figs. 15-18. It can be seen that the largest differences 
occur close to the frontal region for the u field, but hey are almost uniformily 
distributed at the lowest levels in the 0 field for experiment SLSI15. 
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Figure 15. Differences in u component field (m/s) for experiment SLSI23 at t = 36 
hours. 




y(km) 

Figure 16. Differences in u component field (m/s) for experiment SLSI15 at t = 36 
hours. 
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Figure 17. Differences in 0 field (°K) for experiment SLSI23 at t = 36 hours. 




Figure 18. Differences in 0 field (°K) for experiment SLSI15 at t = 36 hours. 
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The second set of comparisons used as reference the same solutions obtained for 
experiment SLHR02 at time t = 36 hours. The frontal width at that time in the high 
resolution experiment was 99 5 km For those SLSI experiments where the time filter 
controlled the high frequency oscillations, the frontal width d- 100 km was uniquely 
related to a certain time. If the solutions had no error we could suppose that for a 
certain d-vaJue the velocity components and temperature fields should be the same as 
those obtained from the control run. Thus, the geophysical fields for the times 
corresponding to a d-value approximately equal to 100 km were chosen for the second 
set of comparisons. Table 7 gives the RMS differences obtained and the corresponding 
frontal widths and times used for each comparison. The values obtained show that the 
frontogenesis process is slower for longer time steps, since the frontal width of 100 km 
is reached at a later time for larger A/. The RMS differences for t/and 0 are again 
relatively small compared to the magnitude of the values obtained in the reference 
solution. The differences for v and w are still relatively large with respect to the 
reference values, as expected. Figures 19-22 show the grid point differences of the u 
and 0 fields for experiments SLSI23 and SLSI 15. 

The figures show that the largest differences occur close to the front. This 
characteristic of semi-Lagrangian schemes in which the large errors are confined around 
the scale collapse region was observed by Kuo and Williams ( 1990). Such a property 
is desirable, since the regions away from the front will not undergo a large effect of 
errors generated close to the region where strong gradients are present. 

Figures 23-29 show the temperature distribution at the lowestcomputational level 
(z = 167 m) for the experiments listed in Table 5 (solid line) for comparison with the 
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values obtained from the control run (dashed line) It can be seen that as the time step 
increases the differences increase but the largest errors remain close to the frontal 
region. 

TABLli 7 RMS DIFFERENCES FOR d - 100 km. 





RMS DIFFERENCES 






Experiment, 
At(sec ), y 


u(nvs) 


v(m/s) 


w(m/s) 

*10 U ’ 


ei°K) 


d(km) 


lime 

(hours) 


SLSI03, 
180, 0 03 


0.667 


0.398 


2.91 


0.215 


101.0 


38.0 


SLSI06, 
360, 0.03 


0.683 


0.422 


3.01 


0.240 


101.4 


38.2 


SLSI09, 
600, 0.03 


0.731 


0.468 


3.24 


0.280 


1008 


38.5 


SLSI12, 

1200,0.05 


0.914 


0.593 


3.78 


0.388 


101.2 


39.3 


SLSI15, 
1800, 0.05 


1.031 


0.622 


3.92 


0.452 


101.5 


39.5 


SLSI19, 
2400, 0.07 


1.155 


0.701 


4.26 


0.531 


101.5 


40.0 


SLSI27, 
3600, 0.10 


1.197 


0815 


5.47 


0.588 


103 


400 
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y(km) 



Figure 19. Differences in u component field (m/s) for experiment SLSI03 at time 
corresponding to d- 100 km. 




Figure 20. Differences in u component field (m/s) for experiment SLSI15 at time 
corresponding to d-100 km. 
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y(km) 



Figure 21. Differences in 8 field (°K) for experiment SLSI03 at time corresponding to 
d- 100 km. 




y(km) 



Figure 22. Differences in 6 field (°K) for experiment SLSI15 at time corresponding to 
d-100 km. 
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For this model the Courant number o is calculated using the following expression: 



0 = (< v Lc +l ^ 1 +,C J ) 

LLULA u mflT « max 



max 



Af 

(Ay/2) 



(5.3) 



where c^^ is the phase speed of the fastest internal gravity wave. An important result 
from this study is that for the SLSI experiments, o ranged from 0.5 (A/ = 180 sec.) to 
10.4 (A t - 3600 sec.) , as opposed to the CFL stabilty criterion that would require a to 
be less than or equal to 1, and the values of the RMS differences obtained for wand 0 
in the experiments with large time steps remained relatively small. This result indicates 
that the gain in computational efficiency does not affect significantly the accuracy of the 
solutions obtained with the SLSI scheme. 




Figure 23. Temperature disturbance at z = 167 m, for experiment SLSI03 (solid line) 
and control run (dashed line), at time corresponding to d- 100 km. 
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Figure 24. Same as Fig. 23 except for experiment SLSI06. 




Figure 25. Same as Fig. 23 except for experiment SLSI09. 
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Figure 26. Same as fig. 23 except for experiment SLSI12. 




Figure 27. Same as Fig. 23 except for experiment SLSI15. 
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Figure 28. Same as Fig. 23 except for experiment SLSI19. 




Figure 29. Same as Fig. 23 except for experiment SLSI27. 
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VI. SUMMARY AND CONCLUSIONS 



In this study a numerical model based on the hydrostatic Boussinesq equations 
was used to simulate atmospheric frontogenesis driven by an irrotational non-divergent 
deformation wind field. The semi-Lagrangian technique was employed to integrate 
numerically the prognostic equations. The model assumes periodic boundary 
conditions in the cross-front direction and the rigid lid assumption is used as the upper 
bounda ry condition . The model also neglects along front variations. A basic sinusoidal 
deformation flow is introduced as the forcing of the frontogenesis process. Experiments 
were performed using the semi-Lagrangian technique associated with two different time 
schemes: explicit and semi-implicit. In the semi-Lagrangian explicit case (SLEX) 
bicubic splines were used to interpolate the variables in space. In the semi-Lagrangian 
semi-implicit case (SLSI) the variables were interpolated in the y-direction only, using 
one-dimensional splines. A frontal width parameter (d-vaJue) was defined at the lowest 
computational level to describe the evolution of the frontogenetical processes. The 
experiments with the SLEX technique were successful in replicating the formation of a 
realistic front in approximately 40 hours. Different time steps were used for the 
integration of the model and solutions which were both numerically and physically 
consistent were obtained with the SLEX scheme for values of Courant numbers as large 
as 1.7. 

The experiments with SLSI scheme presented high frequency oscillations that were 
eliminated by using a time filter. However, no small scale spatial oscillations were 
observed in the solutions. Different filtering effects were tested by changing the value 
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of the filtering parameter y The solutions obtained with the SLEX scheme were not 
affected by the time filter, whereas the SLSI solutions showed that stronger filtering was 
necessary for larger time steps. 

The SLSI scheme showed to be more successful in handling the scale collapse 
process than the SLEX scheme, since the fronts reproduced by the former had 
minimum widths smaller than those obtained by the latter. The SLSI scheme presented 
a tendency of slow ing down the frontogenesis process for increasing time steps 

Experiments were performed with the SLSI scheme with time steps as long as 3600 
seconds, corresponding to a Courant number greater than 10. The solutions obtained 
were perfectly stable, and the accuracy was not significantly degraded, even for longer 
time steps. 

The SLSI solutions also had the characteristic of constraining the largest errors 
close to the scale collapse region. Such a characteristic is desirable since the errors will 
have a smaller impact on the solutions in regions away from the front. 

The results obtained in this study suggest that the SLSI technique is appropriate 
for mesoscale regional models since it is computationally efficient and produces accurate 
results. A different formulation of the scheme where two-dimensional interpolation of 
variables were allowed should be studied to verify if there would be a significant positive 
impact on the accuracy of the solutions. Also, effects of surface topography and 
physical processes like advection of the frontal system, friction, vertical shear of the 
basic wind field and moisture should be investigated in future studies. 
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